'''
Created on Sep 3, 2010

@author: oabalbin
'''
import os
import subprocess
from optparse import OptionParser
from collections import deque, defaultdict
from bx.intervals.intersection import Interval, IntervalTree

class myInterval:
    def __init__(self,chr, start, end, name, score, strand):
        self.chr, self.start, self.end, self.name, self.score, self.strand = \
        chr, int(start), int(end), name, score, strand
        


def get_flow_cell_info(sampleInfoFile,flowCellInfo):
    
    sampleType,tissueType,QC=flowCellInfo[0],flowCellInfo[1],flowCellInfo[2]
    sampleDict = defaultdict(list)    
    for line in sampleInfoFile:
        fields = line.strip('\n').split('\t')
        if sampleType==fields[4] and tissueType==fields[5] and QC==fields[12]:
            sampleDict[fields[0]] = fields[1:]
    
    return sampleDict


def extract_readsfrom_BAM_file(flow_cell_path, analysis_dir_path, path_to_samtools, bed, outputfile):

    # go to flow cell dir
    # This can maybe rewrite using the pysam library.
    # This will allow to make quality control on the fly.
    ### Also Check waht you need to do in order to output the arf and the fastaq file for mirDeep from this point
    ### So you do not need to use the bwa_converter after going again over the files. 
    os.chdir(flow_cell_path)
    samtools_paramters = [path_to_samtools, "view","-h","tophat_aligned_reads.bam",bed]
    subprocess.Popen(samtools_paramters,stdout=outputfile)
    
    # go back to analysis directory
    #os.chdir(analysis_dir_path)


def create_interval_tree(bed_file):
    
    intervals_bed_file = open(bed_file)
    intersecter = IntervalTree()
    
    for line in intervals_bed_file:
        fields = line.strip('\n').split('\t')
        nfields = len(fields)
        if nfields < 3:
            print "The bed file has not enough information"
        elif nfields == 6:
            chr_tmp=fields[0].replace('chr','')
            name, score, strand=fields[3],fields[4],default_strand(fields[5])
            intersecter.insert_interval( Interval( int(fields[1]), int(fields[2]), value={'chr':chr_tmp, 'name':name,'score':score}, strand=strand ) )
        elif nfields < 6:
            chr_tmp=fields[0].replace('chr','')
            name, score, strand="noname","noscore","+"
            intersecter.insert_interval( Interval( int(fields[1]), int(fields[2]), value={'chr':chr_tmp,'name':name, 'score':score}, strand=strand ) )
    
    intervals_bed_file.close()                
    return intersecter

def default_strand(x):
    if not x:
        print x
        return '+'
    else:
        return x


def create_list_of_intervals(intervals_bed_file):
    
    intervals_bed_file = open(intervals_bed_file)
    interval_list=deque()
    for line in intervals_bed_file:
        fields = line.strip('\n').split('\t')
        nfields = len(fields)
        if nfields < 3:
            print "The bed file has not enough information"
        elif nfields == 6:
            chr_tmp=fields[0].replace('chr','')
            name, score, strand=fields[3],fields[4],default_strand(fields[5])
            interval_list.append((myInterval(chr_tmp, int(fields[1]), int(fields[2]), name, score, strand)))
        elif nfields < 6:
            chr_tmp=fields[0].replace('chr','')
            name, score, strand="noname","noscore","+"
            interval_list.append((myInterval(chr_tmp, int(fields[1]), int(fields[2]), name, score, strand)))
    
    intervals_bed_file.close()
    return interval_list

    
def intersect_intervals_tree(list_of_intervals, interval_tree):
    
    intersected_intervals=[]
    for thisintv in list_of_intervals:
        intesected_intervals= interval_tree.find(thisintv.start,thisintv.end)
        # Specific for the repeat program. It could be eliminated from other scripts, AluSg7|Alu|SINE, MIR|MIR|SINE
        repeat_type = thisintv.name.split('|')[1]
        if repeat_type == "MIR": #MIR #"Alu"
            for intv in intesected_intervals:
                if (intv.value['chr']==thisintv.chr) and (intv.strand == thisintv.strand):
                    intersected_intervals.append(intv)
                    #print intv, thisintv.name
    
    return intersected_intervals
    
        
def create_SAM_summary(data_flow_cell_path,analysis_dir_path, path_to_samtools,outputfile_path, myintervals, flow_cell_dictionary):
    
    
    for flow_cell in flow_cell_dictionary:
        outputfile = open(outputfile_path+flow_cell+'_test_tu_complilation.sam','w+a')
        flow_cell_path=data_flow_cell_path+flow_cell+'/'
        print flow_cell_path
                
        for intv in myintervals:
            location = 'chr'+intv.value['chr']+':'+str(intv.start)+'-'+str(intv.end)
            extract_readsfrom_BAM_file(flow_cell_path, analysis_dir_path, path_to_samtools, location, outputfile)
        
        outputfile.close()
    
    


if __name__ == '__main__':
    
    '''
    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-f", "--annotFile", dest="annotFile",
                            help="annotation file for all files to use")
    optionparser.add_option("-d", "--outfolder", dest="outfolder",
                            help="output folder")
    
    (options, args) = optionparser.parse_args()
    '''
    
    sampleInfoFile = open('/data/projects/mirnas/sampleinfo_2010_Aug_08.txt')
    TUs_with_repeats='/data/projects/mirnas/mirnas_discovery/repeat_hits.pad50.bed'
    TUs_aluexplosion='/data/projects/mirnas/mirnas_discovery/2010_09_01_14_30_classifier_transcripts.bed.txt'
    #data_cell_path='/data/projects/mirnas/mirnas_discovery/'
    data_cell_path='/exds/tuxedo/version_003_2010_Mar_26/data/'
    analysis_dir_path='/data/projects/mirnas/mirnas_discovery/'
    path_to_samtools='/exds/sw/bioinfo/alignment/samtools/current/samtools'
    outputfile='/data/projects/mirnas/mirnas_discovery/sam_files_MIR/'
    flowCellInfo = ['Tissue','Prostate','PASS']
   
    #==== main =====#
    flow_cell_dictionary = get_flow_cell_info(sampleInfoFile,flowCellInfo)
    #flow_cell_dictionary={}
    #flow_cell_dictionary['myflowcell']='thiscell'
    trans_tree1 = create_interval_tree(TUs_aluexplosion)
    interval_list = create_list_of_intervals(TUs_with_repeats)
    intersected_intervals = intersect_intervals_tree(interval_list, trans_tree1)
        
    create_SAM_summary(data_cell_path,analysis_dir_path, path_to_samtools,outputfile, intersected_intervals, flow_cell_dictionary)




